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Abstract 

It is shown numerically, both for the two-dimensional Navier-Stokes 
(guiding-center plasma) equations and for two-dimensional magnetohydrodynamics , 
that the long-time asymptotic state in a forced inverse-cascade situation is 
one in which the spectrum is completely dominated by its own fundamental. 

The growth continues until the fundamental is dissipatively limited by its 
own dissipation rate. 
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1, Introduction QP POOR 

It is by now well known that there exist circumstances under which 
the equations of some nonlinear continuous media exhibit turbulent solutions 
which transfer a global quantity to long wavelengths at a rate comparable to 
that at which the quantity is supplied. This is a reversal, for that quantity, 
of the older Kolmogorof f-Richardson picture of spectral transfer to short wave- 
lengths and dissipation there. The first historically important case was the 
case of two-dimensional Navier-Stokes {2D NS) turbulence as initiated by 
Onsager (1949) > Fj?$rtoft (1952), Kraichnan ( 196 T), Batchelor ( 1969 ), Leith (I 968 ), 
and Lilly (I 969 ). The relevance of these calculations for plasma physics 
follows from the fact that the mathematics of the 2D NS case is identical to 
that of the electrostatic guiding-center plasma (e.g., Montgomery 1975). 

The demonstration of similar long-wavelength transfer in three- 
dimensional, incompressible magnetohydrodynamics (3D MHD) was due to Frisch 
and collaborators (Frisch et al 1975 Pouquet et al 1976, Meneguzzi et al 
1981 ). It was soon demonstrated that still a third quantity could be inversely 
transferred to long wavelengths in two-dimensional, incompressible magnetohydro- 
dynaraics (2D MHD) by lyfe et al (1976, 1977a, b) and Pouquet (1978). A review 
of the literature up to I 98 O appears in Kraichnan and Montgomery (I 98 O). 

Because of the difficulty of doing two-dimensional experiments, the 
cited work has all been theoretical and/or computational. The emphasis has 
been on power-law behavior for wavenumber spectra, following the influential 
conjecture of Kraichnan (I967) concerning two possible inertial subranges 
transferring energy and enstrophy in opposite directions. Some investigations 
of dissipation scale behavior have al.so been reported for the 2D MHD case 
(Orszag and Tang 1979, Matthaeus and Montgomery, I 98 O, I 98 I, Matthaeus 1982, 
Frisch et al, I 983 ). 
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onmu.L pfx-: 

OF POOr^ QUALITY 

Pover laws are only obtained persuasively from dimensional analysis. In order 
for them to be seen clearly, they are thought to req.uire very large Reynolds 
numbers (mechanical and/or magnetic) in order that the inertial subranges be 
cleanly separated from each other and from the dissipation range. High spatial 
resolution is req,uired at these large Reynolds numbers, and it is accm’ate to 
say that no reported numerical solution has contained a wide enough range of 
wavenumbers to satisfy the chain of inequalities required for a satisfactory 
test of the power-law behavior in any putatively inverse-cascading situation. 
There is no question, however, that the qualitative effect of large back-trans- 
fer, at least not inconsistent with the proposed power laws, has been seen in 
numerical solutions of the relevant dynamical equations in all three cases 
(Lilly 1969> Fyf^ 19TTb, Meneguzzi et al I 981 ) . 

Theoretical predictions typically have been for unbounded systems, 
so that there were arbitrarily many additional octaves in wavelength for the 
spectrum of the injected and inversely-cascaded quantity to expand into at 
the lower end. Numerical solutions, however, always take place on periodic 
or finite grids, and thereby have associated with them a maximum wavelength 
(minimum wavenumber) beyond which the inversely-cascaded spectrum cannot go. 
Because of the slow time scales involved in the long-wavelength eddy turnover 
times and because of limited computer budgets, computations have not been report- 
ed for times long compared to the time required for the longest allowed wave- 
length to fill up. Both theory and computation have remained vague on the 
question of what happens in an inversely-cascading situation. 

Here, we address this question numerically for both the 2D NS and 
2D MHD cases. In both cases, the fields are driven by a prescribed random 
function, band-limited in wavemmiber space, which serves as the source of the 
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excitutlonH that, are trauai'errcd eiaowhero In the apaoe. The reaoLutien of 
the eomputatien ia net hi(^h aiul the Reynolds number a are not large. There is 
no pretense that the iuertlal oubrangeo are well separated or that power-law 
prediotions are verifiable. It is not noeeasarjs however, to be in thia high 
Reynolds number regime in order to see large amount^i of back transfer, or to 
address the question of wlmt happens to the baek-tr ana f erred exeitatlone. when 
they liavo no lower pXaee in wavenumber spaee to which to migrate. It ia toward 
this latter, qualitative question that the following material ia directed. 

There would seem to be at least three possible scenarios for what 
might oeeur after tiie k spectrum fills in between the forcing baud and the funda- 
mentuti lowest wavenumber (1) The inverse spoetral transfer iiiight cease. 

(a) The rolevaait parts of the k spectrum might all rise together, posaibI.,Y main- 
taining a jxswer law behavior, vltli perhaps a ’’hOvaling region” Just below the 
forcing band, (i) The fumhuiiental k^^j^ might continue to absorb the supplied 
excitations, running off and leaving the rest of the spectnuti, until it is 
limited by its o\fu dissipation rate, These three imw not exhaust the possibili- 
ties, but wliat we wisl\ to demonstrate in the following pages is that alternative 
(3) is in fact Aduit occurs, both for the GD Nf' and dP MUP cases. Wo are unaware 
of miy previous predictions or tests of thia effect. 

In section 2, wc outline tlie computational proeedure by which the 
tests are carried out. In sections i find t, respectively, results are reported 
for tlie dp R'3 and OP MUP cases. Section 5 describes some aemi-qiu\ntitativo 
attempts at modelling iilgebraically the k^^^j.^-domlnated final state, and briefly 
svumnarines ow results. Questions of accui'aey of the mMorlcal method are 
addressed in the Appendix, 
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2. Numerieal Praaedure 

The OD MltP equations can be written in a standard dimensionless 
form which, by the omission of several terms, reduces to the 2D case. 
Specificiilly, 


1^ = -y . Vw + g • VJ + vV^oj + F (1) 

1^ * >y-Va+yV^a*^G (2) 

(see, e.R., I^yfe et al, 1976, iPT7a,b). All fields are Itinctions of x,y mid 
the time t, mid turn indoj^endout of the 2 -coordlnate. The vorticity is 

f') 

CO K -V'lji, wliere ijj is the stream ituiction mid 'tlie velocity y « Vtl/xe^,, The mag- 
netic vector potential is ao„ mid the magnetic field is B = VaxS , The electric 

u "** Zi 

current density has only a a component mid la j = -V^a. The vector fields have 

only X and y components, v == (v„, v , O) mid B = (B„, B , O) , though a constmit 

~ X y ~ y 

dc magnetic field added to B without altering equations (l) 

and (2). Deleting equation (2) and setting B = 0 leaves the 2D N8 system. 

f and G are rmidom forcing f mictions of x, y mid t whose statistics 
will be described below. 

Rectangular jHipiodic boiuidsu’y conditions are assumed, for all fields, 

over a square box of edge 2'ir. All fields are expressed as Fourier series with 

waveiuuubera k = (k^, k^) with integer components. For example, for the vorti- 

city, to = Sj^co(k,t)exp{ik‘ 3 c) . Finite computational x*eaoui'ses make it necessary 

to ti’uncate the Fourier series at the I'elatively low value of tkl = k - l6, 

' ~ ' max max ’ 

though some runs with k = 32 have also been cai'ried out. TVie mechanical 
° max 

. -1 -1 

and magnetic Reynolds numbers, \) and p , must be diosen small enougli to 

suppress wavenumbers significmitly above k , for the results to be physical. 

max 

The V, p chosen for the runs reported were never smaller tlian 0.002. 
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Solution of equations (l) and (2) is by the by-now standard Orszag- 
Patterson spectral method (Orszag 1971 » Patterson and Orszag 1971 » Gottlieb 
and Orszag 1977) > •’Thich has proved extraordinarily useful in turbulence compu- 
tations. As a time-stepping scheme, we use 


Qn‘^l/^(k) = Q^(k) + P( Q^(k) )At/2, 

(3) 

Q“'*’^(k) = Q^(k) + P(Q^'^^^^(k))At. 

(4) 


Here dQ(k,t)/dt = P(Q(k,t)) stands symbolically for the Fourier-transformed 
version of either equation (l) or (2), and Q(k,t) represents all the Fourier 
coefficients of all the fields. The index n indexes the nth time step of 
duration At. 

The forcing functions F and G have Fourier transforms F(k,t), 

G(k,t) which are chosen at each lattice point of the k space, according to the 
following recipe. Let H’^(k) be the real or imaginary part of either F(k) or 
G(k) at time step n. Then (Fyfe et al, 1977b) » 

H“‘*'^(k) = fH^(k) + /l-f^ J^(k). (5) 

J^(k) is a random number taken from a Gaussian random number generator. The 

(• 

"memory fraction" f is chosen between 0 and 1, and for the runs reported, we 

used f = 0.95 fo^ both the 2D NS and 2D MHD runs. F(k,t) and G(k,t) are set 

identically zero outside a "forcing band" in wavenumber space, k ^ 1 ^F,„ay • 

2 2 

In the riins reported, k« = 55 aad k™ =70. 

^min ^max 

In both the 2D NS and 2D MHD runs, we started with an empty spectrum 
and followed the evolution over a considerably longer interval than previously 
(Fyfe 1977b). For the 2D MHD run, the F and G functions were \in correlated, but 
acted in the same band between k“^ = 55 and 70, with statistically equal 
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magnitudes. The "filling up" of the spectrum at the longest wavelengths is 
a slow process, and consumes much computer time; it is an inherent featxire of 
such computations, and apparently the reason these computations have not 
been carried out previously. 

3. 2D NS Case; Inverse Energy Cascade 

Results for the Navier-Stokes run are shown in figures 1-6. 

1 P 

Figure 1 shows the total energy E(t) = 2E(k) - |y(k,t)| as a 

function of time. The two most significant features of the curve are: (i) the 

total energy increases systematically with time, as is expected for an inverse- 
ley cascading g.uantity; (ii) the total energy appears to be approaching asymp- 
totically to a constant value, whereas for an unbounded system in the limit of 
infinite Reynolds numbers, a linear increase of E(t) with t might be expected. 
As will be discussed in section 5» the curve in figure 1 is rather well fit 

by the expression E = E [l-exp(-t/T) ] with E „ = 3*9 and t = 124. 

max max 

Also in figure 1, we show, as sq.uares, the energy contained in the 

op ^ . 

fundamental k = k . =1. It will be apparent that by t = 30, it already 

min 

accounts for a significant fraction of the total. 

Figure 2 shows the total enstrophy R = Sfi (k) = —S, |m(k,t) | versus 

k ~ iS 

time. The enstrophy shoUld be the directly cascading quantity, and reaches a 
value about which it fluctuates while the energy continues to grow. 

2 

Figures 3a and 3b show time averaged modal energy spectra [y(k)| 

2 

versus k . These are averaged over all the values of k corresponding to a 
given k . Figure 3a shows averages over 1000 time steps ending at t = 29.3, 
and figure 3b shows averages over 1000 time steps ending at t = 273.44. 

p 

(Every fourth point is plotted for k‘^>70, to avoid cluttering the graphs.) 
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The crucial reoult of this paper is thought to be illustrated by the fact that 
the fundcunentol has come to contain about ninety per cent of the total energy. 

Though the separation of the spatial scales is far too scant to 
identify inertial subranges, we show in figures 3'i»b least-sq,uares, straight- 
line fits to the spectra above and below the forcing bonds. In figure 3b, the 
exponents above and below are -l.l6 and - 1 . 62 , corresponding to -1.32 and -2.2h 
for omnidirectional spectra, whereas the Kolmogoroff-style dimensional-analysis 
predictions are -1.6T and -3-0» omni-directional. We attach no particular im- 
portance to these numbers; the power-law predictions are inapplicable. 

Figure 4 displays large-scale stream function contours of constant ijj 
at t = 273.44, showing the results of domination of the spectrum by the funda- 
mental . 


Figiu’e 5 illustrates the behavior of a typical intermediate = 50) 
modal energy, fluctuating about a nearly constant value. This should be con- 
trasted with the continued giwfbh of the fxindamental displayed in figure 1. 

Figure 6 exhibits the Carnevale ( 198 I) entropy, S = E toj y(k,t) | ^, as 
a function of time. The function q.uickly maximizes, and then gradually decays. 
The interpretation of the decay might loosely be that as the system becomes 
dominated by the fundamental, the effective degree of randomness decreases. 

(The action of viscosity is typically to decrease entropy.) This entropy has 
its origins in information theory, and descends from the Burg entropy (Burg, 


1975). 


4. The MHD Case 

The results of the 2D MHD run are shown in figures 7-11. The spec- 
trvim is forced both mechanically and magnetically, with the ratio between F 
and G chosen to supply kinetic and magnetic energies at statistically eq.ua! 


rates. 
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The Toreing bond (55^h s70) was chosen to be the same as in the 
Navier-Stokes case, and for the run we exhibit, p = v - 0.01, with a time 
step t a 1/25(5. MID runs are typically noisier than NS ones, and the com- 

p 

puted Kolmogoroff dissipation wave niunbers tend to bo higher. The (32) 
runs we have performed with the above parameters tend to lead to Kolraogoroff 
dissipation wave numbers significantly greater than Even though the 

behavior exhibited for the 2D MHD case closely parallels the 2D NS behavior, 
for (32) , we have wought to enhance our confidence in the reality of the 
effect by also performing some (64) runs. Because (64) runs for times as 
long as those reported in section 3 are prohibitively expeniLiive and time 
consuming, we have resorted to the following trick. We computed a (32) r-vm, 
starting with on empty spectrum, up to time 97 >(56. At this time, we used the 
computed Fourier coefficients as initial data for i 5000 time step run (up 
to t - 117.19) with a (64)^ resolution. (Both the (32)^ and (64)^ runs were 
performed on a Cyber 203: a 64-bit machine that is essentially without round- 
off error.) Even though the Kolmogoroff dissipation wave number was unresolved 
by k /k . - 16 dui’ing the first interval (typically, the dissipation wave 

number was 33) , it was reasonably well resolved during the second interval 

(typically, the dissipation wave nmnber was also about 33) when k /k . was 

max min 

32. No significant changes were observed from one piece of the run to the 
other; the second piece, regarded as accurate, continued the behavior ob- 
served in the first piece, regarded as possibly inaccurate. 

Figure 7 shows the mean square vector potential vs. time, with 
squares illustrating the contribution of the fundamental k^^^ = 1 (again, 

2 

values shown are averages over all values of k corresponding to a given k ) . 
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The ti’ancltiou from a (33)*“ computation to a computation occurs at 

t = 97.66, Figure 8 chows, for the some run, the total energj*' vs. time; the 
points laheled "B” are the magnetic contriliutions, and the points labeled "V" 
are the mechanical (kinetic) energy. Tlie higher values of magnetic energy 
are typical of forced MID runs. Figures 7 and 8 should he compared with 
figures 1 and 2 for the Navier-Stokes case. 

Tlie vector potential spectrum (|A(k)|^ = [B(k)|^/k^ vs. k^) is shown 
in figui'e 9» averaged over 50 values, spaced 100 time steps apart, ending at 
t - 97.66. Figure 10 shows the vector potential spectrum, averaged over 25 
values, spaced 100 time steps apart, ending at 117 • 19* The dominance of the 
spectriuu by the fundamental should bo compared with that of figure 3b for 
the 2D NS case. Figure 11 shows the instantaneous contour plots ("magnetic 
islands") for vector potential at t = 117.19 • 

In summary, the same qualitative behavior, with the role of the 
energy replaced by mean square vector potential, is seen for the 2D MHD 
case as was observed in the 2D NS case. 

5. Algebraic Model 

The indications are, for the runs reported in sections 3 and h, that 
a steady state is being approached in which the fundamental dominates the spec- 
trum, and the transfer to it is balanced by its own dissipation. This has led 
us to try a crude algebraic model of the dynamics of this process, which we 
illustrate for the 2D NS case. 

Once the fundamental has begun i > dominate, most of the new energy goes 
into the fundamental, so that E= E(kj^n)> ^ assume that the 

total enstrophy has become constant, dfJ/dt = 0, and a significant fraction of the 
enstrophy may be lodged in the values of k > k^^^ ( as the Reynolds numbers increase , 

this conclusion should become more accurate). Let Q - + f2(k.), where 12^ = SI2(k). 

k>k ~ 
rain 
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The evolution of the energy and ens trophy are now modelled hy 


-2vR + fg (6) 

£ = 0 = -2\>zA(^) + k|fg (T) 

where is the rate of supply of energy hy F and kj, is a wavenumher at the 
center of the forcing hand. 

fg may he eliminated from equations (6) and (7) to give 

f = -avsi A(fe) 

■ivjt 

or (noting that = l) > 

~ + 2 (1 - -%) E= -2vn^ + ^ Sk^n(k) (8) 

F h.>kuiij<i 


The right hand side of equation (8) can he evaluated in terms of the spectrum 
lyCk)!*^ = C^k”^! for k helow the forcing hand and Cj_k~'^2 above it, with n^^, ng, 
and taken from the computed values in figure 3h. The result is that equa- 
tion ( 8) becomes of the form 


dE . 1 „ _ max 
dt T T 


(9) 


with E = 3.58 and t = 25^ from the calculation. The solution, E = E [1 
max * max 

-exp(-t/x)] fits the shape of figxire 1 rather well, hut the numbers are not 
accurate, for E =3*9 from figure 1, and x = 124. It might he, in the high 

|TIB*X 

Reynolds number limit, that this simple model in which the energy dissipation 
occurs mainly in the fundamental and balances the input, could he sharpened 
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considerably, but at such a low spatial resolution it has not seemed worth- 
while to do GO, 

A similar development can be carried out for the 2D MHD case, with 

the resulting equation for A predicting dA/dt - -A/t + A /t, and a (com- 

in8<x 

parably inaccurate) fit to the computed behavior can be obtained. 

In stimmary, we believe that we have identified the correct quali- 
tative behavior for the long-time state of the randomly forced 2D NS and 2D 
MHD equations in the presence of periodic boundary conditions. Namely, the 
fundamental continues to absorb the back-transferred quantity until it domi- 
nates the spectrum and may be limited only by its own dissipation. 

Acknowledgments . This work was supported in part by the U-.S, Department of 
Energy and in part by National Aeronautics and Space Administration Grant 

nsg-t4i6. 
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Appendix; Questions of Long-Term Acouracy 

An obvious question which has no simple answer is whether there is 
a sense in which the code con be called accui’ate after the many ed(^ turnover 
times represented in these runs. Accuracy checks for such problems are few, 
and mostly consist of Inferences drawn from conservation laws of the ideal 
invariants ("rugged" invariants) which survive the limitation to finite numbers 
of Fourier coefficients. The fact is that the code will not conserve these 
invariants out to the times we are reporting. These long times are dictated 
by the physics of the situation (it takes that long for the fundamental to come 
to dominate) and there are no totally conclusive statements about the accuracy 
that we can make. 

We can, however, derive conside'^able reassurance fx-om the following 
observatior?;,; , Except for the very beginning (t^ 5, say) there are no periods 
of ‘•ABUi in which large qualitative changes occxxr in the behavior of the computed 
quantities. The code is accurate (in the above sense) over this initial time 
interval, and after that, there is only a slow, systematic evolution, always 
in the same direction. The code can be re-started and run accurately over any 
slice of time during the evolution after randomizing the phases of the Fourier 
coefficients; e.g., 6j(k)-^w(k) exp (ic})) at some instant, where ^ is a random 
number between -tt and it. 

Figure 12 shows the results of such a typical restart at t = 58.6 
for the 2D NS run, in a plot of E(t) vs t. The two sets of random numbers F 
and G are the same sets of random numbers. There is no noticeable departure 
of the two solutions from each other, even after 10,000 time steps. Plotted 
in figure 13 are time averaged modal energy spectra for the two runs, over 
1000 time steps ending at t = 9T-66. Again, no noticeable differences appear. 
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The accumulation of iound-off errors may appear as effectively only 
a slightly different set of random numbers provided by F and G. 

If the problem is linearized and the amplification factor is calcu- 
lated for the numerical scheme, an e-folding time turns out to be about 
t = 2000, whereas ve have integrated only cut about t = 275* 

These three considerations lead us to believe that what we have 
observed are properties of the randomly-forced 2D NS and 2D MHD equations, and 
not properties of bhe numerical scheme. 

We believe the most serious limitation on the accuracy to be connected 
with the spatial resolution. Our computed Kolmogoroff dissipation wave number 

gets as large as 27 for the 2D NS case (k = l6) and as large as 33 for the 

ni8<x 

2D MHD case (k , _ = l6) . For this reason, we have run the additional segment 
max ^ 

of the 2D NS run as described in section 3* As noted there, no reversal of 

any of the observed properties was observed when the dissipation wavenumber 

was < k =32. 
max 
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Figure 1, 

Figure 2. 
Figure 3a. 

Figure 3"b. 

Figure 4. 
Figure 5* 

Figure 6. 
Figure T- 

Figure 8. 
Figure 9. 


Figure 10. 


Figure Captions 

Total energy vs. time for forced 2D NS run. Sq.uares are the con- 

2 2 

trihutions of the fundamental, k = k . =1. 

min 

Total enstrophy vs. time for same 2D NS run shovn in figure 1. 

2 

2D NS modal energy spectrum, averaged over all k for given k , 

2 

vs. k . Averages are taken over 1000 time steps ending at t = 29*30. 

o 

Forcing hand is 55-k <T0. 

Same as figure 3a, except time averages are over 1000 time steps 

ending at t = 273.44. Note dominance of the fundamental, k^ = 1. 

Contours (streamlines) = const, for 2D NS run at t = 273*44. 

Temporal history of a typical (k) modal energy for 2D NS riin. 

2 

Value shown has j kj = 50. 

Garnevale entropy vs. time for 2D NS run. 

Total mean square vector potential vs. time for forced 2D MHD run. 

2 

Squares are contributions of fundamental k =1. Run goes from 
(32)^ to (64)^ at t = 97*66. 

Total energy vs. time for 2D MHD forced run shown in figure 7. 

"B" = magnetic energy; "V" = kinetic energy. 

Spectrum )A(k)j^ = |B(k)|^/k for mean square vector potential for 

2D MHD run shown in figures 7 and 8. Average is over all k's 

2 

corresponding to same k and over 50 values 100 time steps apart 
ending at t = 97.66. 

Spectrum |A(k)|^ = |B(k)|^/k^ for forced 2D MHD run at (64)^, 

averaged over 25 values 100 time steps apart, ending at t = II7.I9* 

2 

Note the dominance of the fundamental k =1. 


Figure 11. Contours of constant a (field linen) at t = II 7 .I 9 foi’ forced 
2D MHD run. 

Figure 12. The result of randomization of phases of a 2D NS run. Two runs 

are identical up to t = ^Q,6 (arrow). Then the phases are random- 
ized for one run and the two hecome distinct beyond that time. The 
departures of the circles from the crosses indicate the divergence 
of the two runs thereafter. 

Figure 13. Modal energy spectra for the two runs shown in figure 12. Note 
that the fundamental (k =1) values are essentially identical. 
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